The "zeroth law" of turbulence: Isotropic turbulence simulations revisited 
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The dimensionless kinetic energy dissipation rate C e is estimated from numerical simulations of 
statistically stationary isotropic box turbulence that is slightly compressible. The Taylor microscale 
Reynolds number (Re^) range is 20 < ReA < 220 and the statistical stationarity is achieved with 
a random phase forcing method. The strong ReA dependence of C £ abates when ReA » 100 after 
which C e slowly approaches « 0.5, a value slightly different to previously reported simulations but 
in good agreement with experimental results. If C e is estimated at a specific time step from the time 
series of the quantities involved it is necessary to account for the time lag between energy injection 
and energy dissipation. Also, the resulting value can differ from the ensemble averaged value by up 
to ±30%. This may explain the spread in results from previously published estimates of C s . 



PACS numbers: 47.27.Ak, 47.27.Jv, 47.27.Nz, 47.27.Vf 
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I. INTRODUCTION 



The notion that the mean turbulent kinetic energy dis- 
sipation rate e is finite and independent of viscosity v was 
originally proposed by G. I. Taylor^. Its importance is 
so recognized now that it is commonly referred to as the 
"zeroth law" of turbulence. Its existence was assumed by 
von Karman and Howarth, Loitsianskii and also, signifi- 
cantly, KolmogorovQ in establishing his celebrated simi- 
larity hypotheses for the structure of the inertial range of 
turbulence. Kolmogorov assumed the small scale struc- 
ture of turbulence to be locally isotropic in space and 
locally stationary in time - which implies the equality 
of turbulent kinetic energy injection at the large scales 
with the rate of turbulent kinetic energy dissipation at 
the small scales Q. Although this view should be strictly 
applied only to steady turbulence, the mechanism of the 
dissipation of turbulent kinetic energy can be considered 
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the most fundamental aspect of turbulence not only from 
a theoretical viewpoint but also from a turbulence model- 
ing viewpoint. Indeed, the mechanism that sets the level 
of turbulent dissipation in flows that are unsteady is a 
difficult, if not intractable, aspect of turbulence model- 
ing. 

The rate of turbulent kinetic energy dissipation is de- 
termined by the rate of energy passed from the large-scale 
eddies to the next smaller scale eddies via a forward cas- 
cade until the energy is eventually dissipated by viscosity. 
Thus, C E defined as, 



C £ =eL/u h \ 



(1) 



(here, e is the mean energy dissipation rate per unit 
mass, L and u' are characteristic large length and ve- 
locity scales respectively) should be independent of the 
Reynolds number and of order unity. An increase in 
Reynolds number should only result in an increase in the 
typical wave number where dissipation takes placeQ. In 
the past few years there have been a number of numerical 
(see Ref. [a and references therein) and experimental (see 
Refs. 0, UT3 f° r recent results) efforts to determine the 
value of C e and its dependence on the Reynolds number. 
Perhaps the most convincing of these are the numerical 
attempts since there is no re-course to one-dimensional 
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surrogacy as there is for experiments. Notwithstanding 
this fact, there is good agreement, both numerically and 
experimentally, with the long held view that C e is ~ 0(1) 
when the Reynolds number is sufficiently high. The col- 
lection of isotropic simulation results for C £ shown in Rcf. 
H indicates that "high enough" Reynolds number ap- 
pears to be Re A ~ 0(100). Here, Re A (= ^[Ib/ve] 1 / 2 ) is 
the Taylor microscale Reynolds number. At higher Re a, 
slow Re a dependencies for C e , such as that proposed by 
LohseQ, cannot be ruled out. Measuring such Re a de- 
pendencies, either numerically or experimentally, will be 
close to impossible. 

One unresolved issue is that raised by Sreenivasan|l(Tj| . 
After assembling all the then known experimental de- 
caying grid turbulence data|ll| and numerical data for 
both decaying and stationary isotropic turbulence he con- 
cludes that "the asymptotic value (of C e ) might depend 
on the nature of large-scale forcing, or, perhaps, on the 
structure of the large scale." He also demonstrates [l2j 
in homogeneously sheared flows that the large structure 
does influence C e . However, it might be argued that these 
results were obtained at low Reynolds numbers and the 
issue of a universal asymptotic value for C e could still 
be considered open. Alternatively it could be argued 
that homogeneous shear flows and the like are strictly un- 
steady turbulent flows and the zeroth law, in its simplest 
guise, should not be expected to apply to such flows. The 
possibility of some characteristics of large-scale turbu- 
lence being universal should not be ruled out. The recent 
observation that input power fluctuations, when properly 
re-scaled, appear universal^] may be construed to sug- 
gest the possibility of universality for C e . The aim of the 
present work is to estimate C £ from direct numerical sim- 
ulations (DNS) of statistically stationary isotropic turbu- 
lence and compare with previously reported DNS results 
(summarized in Fig. 3 of Ref. 0) and experiments car- 
ried out in regions of low (dll/dy ~ dU/dy\ max /2) or zero 
mean shear. The present DNS scheme differs from meth- 
ods already reported in that a high-order finite difference 
method is used. To our knowledge, these are the first fi- 
nite difference results for C e . Hence, it is worthwhile to 
test if different numerics and forcing at the large scales 
result in vastly different values for C e to those already 
reported. 



II. NUMERICAL METHODS 



where Sy = \(uij + Ujj) — hSij V • u is the traceless rate 
of strain tensor. In the numerical simulations the system 
is forced (stirred) using random transversal waves given 
by 



f(x, t) = f e cos [ifc(t) • x + i<j)(t)) , 



(5) 



where k{t) is a wave number with magnitude between 1 
and 2, while <p(t) is a phase between — n and it. Both 
cf>(t) and k(t) are chosen randomly at each time step giv- 
ing a forcing that is delta-correlated in time. The ran- 
dom unit vector e is perpendicular to k and the forcing 
amplitude fo is chosen such that the root mean square 
Mach number for all runs is between 0.13 and 0.15 which 
is not too dissimilar to that found in the wind-tunnel 
experiments to be discussed in the next section. For 
these weakly compressible simulations, the energies of 
solenoidal and potential components of the flow have a 
ratio E pot /E BO i w 10~ 4 -10~ 2 for most scales; only to- 
wards the Nyquist frequency (henceforth fc max ) does the 
ratio increase to about 0.1. It is thus reasonable to as- 
sume that compressibility is irrelevant for the results pre- 
sented here whilst at the same time the present results 
can be considered more comparable and relevant to ex- 
perimental wind tunnel flows than the perfectly incom- 
pressible simulations published so far. The code has been 
validated in previous turbulence studies and the reader 
is referred to Refs. flll ITU and the code web-site[T^ 
for more information. 

The simulations are carried out in periodic boxes with 
resolutions in the range of 32 3 — 512 3 grid points. The 
box size is L x = L y = L z = 2tt, which discretizes the 
wave numbers in units of 1. The viscosity v is chosen 
such that the maximum resolved wave number /c max is 
always greater than 1.5/r], where 77 = (i/ 3 /e) 1 / 4 is the 
Kolmogorov length scale. 

To be consistent with previously published DNS stud- 
ies, the total kinetic energy E is defined as, 

E tot = \ (u 2 ) = \u' 2 = f E(k)dk, (6) 
2 2 J Q 

the integral length scale L is defined, 

(■fcrnax 



L 



2u' z Jo 



k~ x E{k)dk, 



(7) 



and the average turbulent energy dissipation rate is de- 
fined as 



The data used for estimating C e are obtained by solv- 
ing the Navier Stokes equations for an isothermal fluid 
with a constant kinematic viscosity v and a constant 
sound speed c s . The governing equations are given by 



(d t + u ■ V) u = -c 2 V Lap + / visc + / 
(9t +u- V)lnp = -V • u. 



The viscous force is 



fv 



v (V 2 it + § VV ■ u + 2vS ■ Vlnp) 



(2) 
(3) 

(4) 



2v 



k 2 E(k)dk. 



(8) 



Angular brackets denote averaging over the box volume. 
After each run has become statistically stationary (typi- 
cally 1-2 eddy turnovers T = L /u') the average statistics 
are estimated for the remaining total run time. Table ^ 
summarizes the average statistics for each run. Compar- 
ing Runs C and D in Table |l] indicates that there is little 
difference in the average C £ for simulations resolved up 
to 77fc max = 1.5 from r/fc max = 3. 
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Run 


N 


Re a 


Ttot/T 


!^(X10 4 ) 


e(xl0 5 ) 


At/t K (xlO 2 ) 


L 


A 


u'(xW 2 ) 


''"max / T 




?7(xl0 3 ) 


^max 


A 


32 


20 


31 


40 


24 


1.9 


1.9 


1.2 


7.1 


0.15 


1.2 


128 


2.1 


B 


64 


42 


30 


15 


22 


1.5 


1.6 


0.81 


7.8 


0.37 


0.75 


63 


2.0 


C 


128 


90 


11 


4.0 


24 


1.5 


1.3 


0.43 


8.4 


0.62 


0.54 


23 


1.5 


D 


256 


92 


19 


4.0 


21 


0.71 


1.4 


0.45 


8.1 


0.69 


0.53 


24 


3.0 


E 


256 


152 


20 


1.6 


21 


1.1 


1.4 


0.29 


8.4 


0.74 


0.49 


12 


1.5 


F 


512 


219 


7 


0.80 


25 


0.86 


1.3 


0.20 


8.9 


0.67 


0.47 


7 


1.7 



TABLE I: Examples of DNS parameters and average turbulence characteristics. N is the number of grid points in each of the 
Cartesian directions, Re^ is the Taylor microscale Reynolds number = u'X/v, T tot is the total run time after the run became 
statistically stationary, T is the eddy turnover time = Lju , At is the run time increment, t K is the Kolmogorov time scale 
= ^ 1/ ^ 2 e _1//2 , A is the Taylor microscale = u' 'y/ lhv/s, r max is the average time for the energy cascade from large to small scales, 
and r\ is the Kolmogorov length scale = ^ 3 ' /4 £ _1/ ' 4 . 



III. RESULTS 



Numerical results 



In this section results for the higher order finite dif- 
ference numerical simulations are presented. The simu- 
lations began with N = 32 3 and each subsequent larger 
box size began with a velocity field interpolated from the 
previous box size. Figures^a)-(d) show example time se- 
ries from Run E (N — 256 3 ) for the fluctuating velocity 
u, the fluctuating integral length scale L, the fluctuat- 
ing kinetic energy dissipation rate e and the fluctuating 
Reynolds number ReA respectively. Initially, the turbu- 
lence takes a short amount of time to reach a statistically 
stationary state - a consequence of stabilizing the new 
run from the previously converged run. The fluctuating 
quantities shown in Figures ^a)-(d) are not unlike those 
encountered in a wind tunnel. Indeed, Fig. da) could 
easily be mistaken for a hot-wire trace of a turbulent 
flow. This is in stark contrast to some pseudo-spectral 
methods that use negative viscosity to maintain a con- 
stant energy level. For example, it is worth comparing 
Figs, nfb)-(d) with those shown in Ref. .18] [i.e. their 
Figs. (2)-(7)]. The pseudo spectral results show that the 
same quantities only fluctuate with a comparatively long 
period. 

Given that the statistics are fluctuating, although they 
are statistically stationary, it is tempting to plot the in- 
stantaneous C E as a function of ReA ■ Figure [21 shows C e 
calculated in such a way. The ReA dependent trends are 
obviously not as expected. However, it is worth noting 
the apparent range for C e when ReA > 50 is « 0.3 — 0.7 
which is the range of previously published DNS results. 
This may explain the scatter in previously published DNS 
results if C e is calculated from a subjective choice of e, L 
and u' at a single time step e.g. as in Ref. |||. The reason 
for the incorrect ReA dependence for C £ can be gleaned 
from Figs. [Ha) and (b). Figure [Ha) shows that an in- 
tense burst in turbulent kinetic energy u 2 (an example 
is noted by the arrow) can be observed some maximum 
time lag r max later in the turbulent kinetic energy dissi- 
pation rate [Figure ^b), again noted by an arrow]. More 
about the significance of r max will be discussed later in 



Section llll HI By noting that there is a strong correlation 
between intense events of u 2 and L on the one hand and 
e on the other hand it is possible to estimate r max from 
the maximum in the correlation between u' 3 / L and e by 



Pu' 3 /L,e( T ) 



[u>%t)/L(t)][e(t + T)] 



u' 3 (t)/L(t) e{t 



(9) 



Figure [21 shows an example for Run E. The maximum 
time lag T max corresponding to the maximum in p u ^y L e 
is indicated by the up arrow j. 




16 20 



FIG. 1: Example time series from Run E, N = 256 3 , average 
ReA ~ 152. (a), u'\ (b), s; (c), L; (d), ReA. Here, the eddy 
turnover time T — L/u . The up arrows f indicate correlated 
bursts of v! and e. 



4 



With this done for all runs it is possible to shift the 
time series of e(t) for each run by its respective r max 
and correctly calculate the instantaneous magnitude of 
C e . Figure [5] shows the newly calculated ReA dependence 
of C e using the correct time lag r max for each of the 
runs. A number of comments can be made. Firstly, the 
dimensionless dissipation rate C e appears to asymptote 
when ReA > 100. The asymptotic magnitude C e w 0.5 is 
in good agreement with the consensus DNS results pub- 
lished so far i.e. C e » 0.4 to 0.5. (see Ref. and ref- 
erences therein). Having said this and given the present 
demonstration that it is incorrect to estimate C £ from 
a single time snap shot it would be interesting to recal- 
culate previously published results based on subjective 
choices of the quantities involved for estimating C e by 
using the entire time series. Lastly, the present results 
verify the use of a high-order finite difference scheme and 
also prove that the zeroth law applies to slightly com- 
pressible turbulence. 




Re, 

FIG. 2: Incorrectly estimated C e as a function of ReA- +, 
Run A; V, Run B; x, Run C; □, Run D; o, Run E; A, Run 
F. Ensemble averages can be found in (Table [J. 



Having estimated r max and assuming it approximates 
the average time r for the energy to cascade from the 
large energetic scales to the small dissipative scales it is 
worth comparing the present results with a simple cas- 
cade model such as that discussed by Lumley^]. Using a 
forward cascade model, whereby the spectrum is divided 
logarithmically into eddies which have the same width in 
wave number space as their center wave number, the to- 
tal time taken for energy to cross the spectrum, assuming 
that all energy is passed directly to the next wave num- 
ber, 

t = r max = 2(~) (l - 1.29^/l5/[Re£C e ]). 
Here, we have substituted (15/[Re^C e ])s for Lumley's 

— 1/2 

large scale Reynolds number dependence Re L . In non- 




FIG. 3: An example of the correlation p u /3/ Le , Eq. for 
Run EN — 256 3 . The up arrow f indicates the location of 
W/T « 0.74. 



dimensional form, 

r+ =2(l-1.29 v /l5/[R*£C e ]). (10) 

As noted by Lumley, little attention should be paid to 
the numerical values of the coefficients, though attention 
should be paid to the exponent for ReA- For small val- 
ues of t + , e.g. t+ < 1, the large scale energy is directly 
affected by viscosity and has little chance of transferring 
energy in a classical cascade manner, whilst for large val- 
ues of r + , e.g. t + > 1, the large scales have the time 
for grinding down energy unaffected by viscosity. The 
asymptotic assumption of Eq. (|1(J|) is 2 eddy turnovers. 
Figure Wi a ) shows the ReA dependence of r + compared 
with Eq. I|10|) . The present results are much lower than 
the prediction of Eq. (|10|l and this is probably indicative 
of the fact that the energy cascade is not a simple full 
transfer of energy between neighboring wave numbers, for 
low ReA at least. It is more likely that, whilst most of the 
energy is passed to neighboring wave numbers, a dimin- 
ishing amount of the energy is passed to all higher wave 
numbers. What is noticeable from the present results is 
that t + 1=3 1 will not occur until Rca ~ 300 which is an 
ReA at which the ReA dependence of C e will become, ei- 
ther numerically or experimentally, unmeasurable. There 
is no reason not to expect that at high enough ReA full en- 
ergy transfer may occur between neighboring wave num- 
bers. Using Eq. JTDJ , Fig. HO indicates that not until 
Re A - 0(1O 3 ) will t+ « 2. 

B. Experimental results revisited 

Results for the present experiment, originally pub- 
lished in Ref. |y] , are updated here with more data within 
the range 170 < ReA < 1210. Detailed experimental con- 
ditions can be found in Refs. 0, Q and need not be re- 
peated here. The main group of measurements are from a 
geometry called a Norman grid which generates a decay- 
ing wake flow. The geometry is composed of a perforated 
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plate superimposed over a bi-plane grid of square rods. 
The flow cannot be classed as freely decaying as the ex- 
tent of the wind tunnel cross section (1.8 x 2.7 m 2 ) is 
approximately 7 x 11 L 2 . For all the flows presented in 
Ref . |6j , signals of the fluctuating longitudinal velocity u 
are acquired, for the most part, on the mean shear profile 
centerline. For the Norman grid, some data is also ob- 
tained slightly off the center-line at a transverse distance 
of one mesh height where dll/dy ~ dU /dy\ meix /2. 
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FIG. 4: Example of the offset time series for Run E (t„ 



150. Note that the peak 
-, u' 3 /L(t/T)---- e ([t- 



0.74), N = 256 . , average Re a 
events are now well correlated. - 

Tmaxl/T). 



All data are acquired using the constant temperature 
anemometry (CTA) hot-wire technique with a single- 
wire probe made of 1.27/zm diameter Wollaston (Pt-10% 
Rh) wire. The instantaneous bridge voltage is buck-and- 
gained and the amplified signals are low-pass filtered fi p 
with the sampling frequency f s always at least twice 
fi p . The resulting signal is recorded with 12- bit reso- 
lution. Time lags r and frequencies / are converted to 
streamwise distance (= tU) and one-dimensional longi- 
tudinal wave number k\ (= 2wf /U) respectively using 
Taylor's hypothesis. The mean dissipation rate e is es- 
timated assuming isotropy of the velocity derivatives i.e. 
e = ei SO = 15v ((du/dx) 2 ). We estimate ((du/dx) 2 ) from 
the average value of Eio(kx) [the 1-dimcnsional energy 
spectrum of u such that u 2 — J °° EiY)(k\)dki and from 
finite differences {(du/dx) 2 ) = (u l+ i - u/) 2 / (U f s ) 2 ]. 

No corrections for the decrease in wire resolution asso- 
ciated with an increase in Re> are made since all meth- 
ods known to us rely on an assumed distribution for 
the three-dimensional energy spectrum. For most of the 
data, the worst wire resolution is « 2r] where r\ is the 
dissipative length scale = v i ^e^J A . The present in- 
vestigation is limited to one-dimensional measurements 
and suitable surrogates for Eq. Although caution 

should be exercised when higher-order moments of a 
one-dimensional surrogate are substituted for the three- 
dimensional equivalent, the use of the mean quantity £i so 
for e should not be too problematic here. The charac- 
teristic length-scale of the large-scale motions L is L p 
and is estimated from the wave number k\ „ at which 




100 150 200 250 300 
Rex 



FIG. 5: Correctly estimated C £ as a function of ReA. +, Run 
A; V, Run B; x, Run C; □, Run D; o, Run E; A, Run F. 
Ensemble averages can be found in (Table Pi. 



a peak in the compensated spectrum kiEm(ki) occurs 



i.e. L p = l/fci )P |J, Il9j. As well the Norman grid 
data, the recent cryogenic decaying grid turbulence of 
White measured using the particle image velocimctry 
(PIV) technique are included. 




1000 



FIG. 6: ReA dependence of inertial range quantities, o, the 
non-dimensional time lag r^ax = r max /T; , Eq. 1101 . 



Figure shows C e for the present data and that of 
Ref. Q. For all of the data, a value of C e ~ 0.5 appears to 
be the average value. It should be noted that estimates of 
C e from the cryogenic decaying grid data are based on the 
transverse equivalents of the quantities that constitute 
Eq. Q). The majority of the scatter for the cryogenic 
data is due to the uncertainty of L which is extremely 
difficult to estimate from PIV data. Figure [7| confirms 
that C e , albeit a one-dimensional surrogate, measured in 
a number of different flows is independent of ReA- It could 
be argued that the rate of approach to an asymptotic 
value depends on the flow e.g. proximity to initial and 
boundary conditions. The asymptotic value C £ « 0.5 
is in excellent agreement with the present DNS results. 
These experimental results are encouraging considering 
that wind-tunnel turbulence is always relatively young 
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compared to DNS turbulence, e.g. the Norman grid 
turbulence has only of the order of 6 eddy turnover times 
in development by the time it reaches the measurement 
station. 




200 400 600 800 1000 1200 
Re, 



FIG. 7: Normalized dissipation rate C e for different experi- 

*8; o, golf ball, 
normal plate, 



mental flows. □, circular disk, 154 < Re A < 1 
70 < Re A < 146; V, pipe, 70 < Re A < 178; 0, 



79 < Re A < 335; A, Norman grid Nl, 152 < Re A < 506; x, 
Norman grid N2 (slight mean shear, dU/dy ^ dU/dy\ ma x/2), 
607 < Re A < 1215, t>, Norman grid N2 (zero mean shear), 
388 < Re A < 1120; <, decaying cryogenic grid turbulence, 
127 < Re A < 3760. 



turbulence and the statistical stationarity is achieved 
with a random phase forcing applied at low wave num- 
bers. The main result of the numerical simulations is the 
demonstration that C e should only be estimated with 
ensemble averaged quantities from the entire time series 
for which the statistics are stationary. If C e is to be 
estimated at each time snap shot it is necessary to cor- 
rectly account for the time lag that occurs from the large 
scale energy injection to the fine scale energy dissipation. 
Even after correctly correlating the energy injection with 
the energy dissipation, the instantaneous value of C £ can 
vary quite considerably (e.g. ±30%) over the extent of the 
simulation. Such a variation may account for the scatter 
in magnitude of C e in previously published results. Both 
the present numerical and experimental results suggest 
that the asymptotic value for C e is « 0.5. In light of this, 
the previously held view that the asymptotic value of 
C e may be dependent on the large scale energy injection 
could be suspect. Lastly, the results presented are strictly 
applicable only to isotropic turbulence that is stationary 
in time. However, it would be interesting to estimate C e 
for simulations of turbulence unsteady in space and/or 
time e.g. anisotropic turbulence or anisotropic homoge- 
neous turbulence with a mean shear because there is lit- 
tle known for these flows about how the turbulent kinetic 
energy is dissipated. 



IV. FINAL REMARKS AND CONCLUSIONS 
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